Powtórzenie

Model GARCH wykorzystywany jest przy modelowaniu i przewidywaniu zmienności w szeregach czasowych. Stosuje się go głównie w analizie finansowej.

Wzór ogólny

\[ GARCH(p, q): \sigma_t^2=\omega+\sum_{i=1}^{p}{\alpha_i\epsilon_{t-i}^2} + \sum_{j=1}^{q}{\beta_j\sigma_{t-j}^2} \]

Najczęściej omawiany przypadek

\[ GARCH(1, 1): \sigma_t^2=\omega+{\alpha\epsilon_{t-1}^2} + {\beta\sigma_{t-1}^2} \]

\(\omega\) - bazowy poziom wariancji warunkowej

\(\alpha\) – współczynnik reakcji na nowe informacje (efekt ARCH)

\(\beta\) – współczynnik pamięci zmienności (efekt GARCH)

\(\epsilon_{t-1}^2\) – kwadrat błędu z kroku \(t-1\)

\(\sigma_{t-1}^2\) – wariancja z kroku \(t-1\)

#| echo: false
#| warning: false
#| error: false
#| output: false
# library(reticulate)
# options(reticulate.python = "/opt/homebrew/bin/python3")
df = pd.read_csv('aapl.us.txt')['Close']

# Liczenie zwrotów procentowych
returns = df.pct_change().dropna()

normal_gm = arch_model(returns, p = 1, q = 1, 
mean = 'constant', 
vol = 'GARCH', 
dist = 'normal')

gm_result = normal_gm.fit()
Iteration:      1,   Func. Count:      6,   Neg. LLF: 13065814402.395237
Iteration:      2,   Func. Count:     19,   Neg. LLF: 194284647014.4725
Iteration:      3,   Func. Count:     32,   Neg. LLF: 2.4201497917599846e+17
Iteration:      4,   Func. Count:     45,   Neg. LLF: 4.565034343361491e+22
Iteration:      5,   Func. Count:     60,   Neg. LLF: 607823.7338190122
Iteration:      6,   Func. Count:     72,   Neg. LLF: 1.1893203852842213e+22
Iteration:      7,   Func. Count:     86,   Neg. LLF: -18810.327116875804
Optimization terminated successfully    (Exit mode 0)
            Current function value: -18810.327071406402
            Iterations: 11
            Function evaluations: 86
            Gradient evaluations: 7

Założenia dotyczące rozkładu standaryzowanych reszt

\[ \text{reszta} = \epsilon_t = \text{predyktowany zwrot} - \text{średni zwrot} \]

\[ \text{standaryzowana reszta} = \frac{\epsilon_t}{\sigma_t} \]

gm_resid = gm_result.resid
gm_std = gm_result.conditional_volatility

gm_std_resid = gm_resid / gm_std


  • Rozkład normalny (domyślna opcja): "normal"

  • Rozkład t-Studenta – grube ogony rozkładu: "t"

  • Skośny rozkład t-Studenta – grube ogony i skośność: "skewt"

Implementacja

normal_gm = arch_model(
  returns, p = 1, q = 1, 
  mean = 'constant', 
  vol = 'GARCH', 
  dist = 'normal'  # założenie o rozkładzie
  )
normal_resid = np.random.normal(0, 1, len(gm_std_resid))

# tstudent
t_gm = arch_model(returns, p = 1, q = 1, 
mean = 'constant', 
vol = 'GARCH', 
dist = 't')

gm_t_result = t_gm.fit()

gm_t_resid = gm_t_result.resid
gm_t_std = gm_t_result.conditional_volatility

gm_t_std_resid = gm_t_resid /gm_t_std

df_hat = gm_t_result.params['nu']
t_resid = np.random.standard_t(df=df_hat, size=len(gm_t_std_resid))
t_resid_scaled = t_resid / np.sqrt(df_hat / (df_hat - 2))

# skośny tstudent
skewt_gm = arch_model(returns, p = 1, q = 1, 
mean = 'constant', 
vol = 'GARCH', 
dist = 'skewt')

gm_skewt_result = skewt_gm.fit()

gm_skewt_resid = gm_skewt_result.resid
gm_skewt_std = gm_skewt_result.conditional_volatility

gm_skewt_std_resid = gm_skewt_resid /gm_skewt_std

df_hat_skewt = gm_skewt_result.params['eta']
skew_hat = gm_skewt_result.params['lambda']

skewt = SkewStudent()
f_skew_resid = skewt.simulate((df_hat_skewt, skew_hat))
skewt_resid = f_skew_resid(len(gm_skewt_std_resid))
Iteration:      1,   Func. Count:      7,   Neg. LLF: 860917.1833719811
Iteration:      2,   Func. Count:     21,   Neg. LLF: 574442.2359680905
Optimization terminated successfully    (Exit mode 0)
            Current function value: -19332.516921943716
            Iterations: 2
            Function evaluations: 28
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      8,   Neg. LLF: 1155974.1387591776
Iteration:      2,   Func. Count:     23,   Neg. LLF: 1440036.9023137486
Iteration:      3,   Func. Count:     36,   Neg. LLF: 758150.776959892
Iteration:      4,   Func. Count:     51,   Neg. LLF: 772107.2252691337
Iteration:      5,   Func. Count:     64,   Neg. LLF: 74349.21108356393
Iteration:      6,   Func. Count:     78,   Neg. LLF: 1114088.0386486105
Iteration:      7,   Func. Count:     93,   Neg. LLF: 1058040.1187104988
Iteration:      8,   Func. Count:    104,   Neg. LLF: 408148.1463894737
Iteration:      9,   Func. Count:    116,   Neg. LLF: 85041.49790235888
Iteration:     10,   Func. Count:    129,   Neg. LLF: 11750.191503565975
Iteration:     11,   Func. Count:    142,   Neg. LLF: -14550.516251339994
Iteration:     12,   Func. Count:    151,   Neg. LLF: 995385.0829001892
Iteration:     13,   Func. Count:    163,   Neg. LLF: 688662.7914191347
Iteration:     14,   Func. Count:    177,   Neg. LLF: 298567.47995581996
Iteration:     15,   Func. Count:    187,   Neg. LLF: 60096.02476851224
Iteration:     16,   Func. Count:    200,   Neg. LLF: 23720.840565413015
Iteration:     17,   Func. Count:    212,   Neg. LLF: 126857.21970662879
Iteration:     18,   Func. Count:    219,   Neg. LLF: -19340.015592720396
Optimization terminated successfully    (Exit mode 0)
            Current function value: -19340.015595069246
            Iterations: 22
            Function evaluations: 219
            Gradient evaluations: 18

Rozkład normalny standaryzowanych reszt

plt.hist(gm_std_resid, bins = 50, 
         facecolor = 'orange', label = 'Standardized residuals')
plt.hist(normal_resid, bins = 50, 
         facecolor = 'tomato', label = 'Normal residuals', alpha=0.7)
plt.legend(loc = 'upper left')
plt.show()

Rozkład t–Studenta standaryzowanych reszt

plt.hist(gm_t_std_resid, bins = 50, 
         facecolor = 'orange', label = 'Standardized residuals')
plt.hist(t_resid_scaled, bins = 50, 
         facecolor = 'tomato', label = 't-Student residuals', alpha=0.7)
plt.legend(loc = 'upper left')
plt.show()

Skośny rozkład t–Studenta standaryzowanych reszt

plt.hist(gm_skewt_std_resid, bins = 50, 
         facecolor = 'orange', label = 'Standardized residuals')
plt.hist(skewt_resid, bins = 50, 
         facecolor = 'tomato', label = 'skewt residuals', alpha=0.7)
plt.legend(loc = 'upper left')
plt.show()

Założenia o średniej

  • Stała średnia (domyślna opcja): "constant"

  • Zerowa średnia: "zero"

  • Średnia modelowana jako proces AR: "AR" (np. AR(1), AR(2), …)

Ocena wydajności modelu

# library(reticulate)
# use_condaenv("/Users/tomasz/opt/anaconda3/bin/python", required = TRUE)

Testowanie istotności parametrów

Test t-Studenta

Dla każdego parametru \(\theta_i\) modelu GARCH przeprowadzamy test istotności, aby sprawdzić, czy jego wartość jest statystycznie różna od zera. W tym celu używa się testu t-Studenta, który sprawdza hipotezę zerową \[ H_0: \theta_i = 0, \] przeciwko \[ H_1: \theta_i \neq 0. \] Statystyką testową jest \[ t_i = \frac{\hat{\theta_i}}{SE(\hat{\theta_i})}, \] gdzie \(\hat{\theta_i}\) to oszacowana wartość parametru, a \(SE(\hat{\theta_i})\) to jego błąd standardowy. Asymptotycznie \(t_i\) ma rozkład normalny, więc wykonujemy sprawdzenie \[ |t_i| > z_{1-\alpha/2}, \] gdzie \(z_{1-\alpha/2}\) to kwantyl rozkładu normalnego standardowego dla poziomu istotności \(\alpha\).

Alternatywnie, możemy patrzeć na poziom krytyczny \(p\)-value, który jest obliczany jako \[ p = 2 \cdot (1 - \Phi(|t_i|)), \] gdzie \(\Phi\) to funkcja dystrybuanty rozkładu normalnego. Jeśli \(p < \alpha\), to odrzucamy hipotezę zerową i uznajemy, że parametr jest istotny statystycznie.

  • W praktyce, zazwyczaj, gdy \(|t_i| > 2\), to możemy uznać, że parametr jest istotny na poziomie \(\alpha = 0.05\).

Testowanie istotności parametrów

Istotność parametrów dla modelu GARCH(1,1)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from arch import arch_model
  • Wczytujemy dane i obliczamy logarytmiczne stopy zwrotu z akcji Google. Wykres logarytmicznych stóp zwrotu przedstawia zmienność cen akcji w czasie.
googl = pd.read_csv('aapl.us.txt', index_col=0, skiprows=2)
googl['log_return'] = np.log1p(googl.iloc[:,0].pct_change())
googl = googl.log_return
googl.index = pd.to_datetime(googl.index, format='%Y-%m-%d')
googl.dropna(inplace=True)

googl.plot(title='Log Returns of Google Stock', figsize=(14, 4))

  • Dopasowujemy model.
model = arch_model(googl, p = 1, q = 1, mean = 'constant', vol = 'GARCH', dist = 't')
model_fit = model.fit(disp='off', cov_type = 'robust')

model_fit.summary()
Constant Mean - GARCH Model Results
Dep. Variable: log_return R-squared: 0.000
Mean Model: Constant Mean Adj. R-squared: 0.000
Vol Model: GARCH Log-Likelihood: 11665.0
Distribution: Standardized Student's t AIC: -23320.0
Method: Maximum Likelihood BIC: -23284.8
No. Observations: 8361
Date: Sat, May 03 2025 Df Residuals: 8360
Time: 15:49:14 Df Model: 1
Mean Model
coef std err t P>|t| 95.0% Conf. Int.
mu -5.9876e-03 3.137e-04 -19.088 3.187e-81 [-6.602e-03,-5.373e-03]
Volatility Model
coef std err t P>|t| 95.0% Conf. Int.
omega 8.3866e-03 5.463e-05 153.514 0.000 [8.280e-03,8.494e-03]
alpha[1] 0.8100 8.104e-02 9.996 1.594e-23 [ 0.651, 0.969]
beta[1] 9.7975e-06 5.900e-03 1.661e-03 0.999 [-1.155e-02,1.157e-02]
Distribution
coef std err t P>|t| 95.0% Conf. Int.
nu 125.9364 1.673 75.284 0.000 [1.227e+02,1.292e+02]


Covariance estimator: robust

Testowanie istotności parametrów

Estymacja parametrów

Maximum Likelihood Estimation (MLE)
Robust Standard Errors (Bollerslev-Wooldridge)

Walidacja założeń modelu

Wykres residuów

  • Szukamy wzorców w resztach modelu, które mogą sugerować, że model nie jest odpowiedni. Wykresy residuów powinny być losowe i nie wykazywać żadnych wyraźnych wzorców.
  • Zgodnie z założeniami reszty powinny być z rozkładu normalnego, widać jednak często występujące ekstremalne wartości.
fig, ax = plt.subplots(2, 1, figsize=(14, 6))
ax[0].plot(model_fit.std_resid)
ax[0].set_title('Standardized Residuals of GARCH(1,1) Model')
ax[1].hist(model_fit.std_resid, bins=30, density=True)
x = np.linspace(-4, 4, 100)
ax[1].plot(x, (1/np.sqrt(2*np.pi)) * np.exp(-0.5*x**2), 'r', lw=2, label='N(0,1)')
#plot t density
from scipy.stats import t
ax[1].plot(x, t.pdf(x, df=model_fit.params['nu']), 'g', lw=2, label='t-distribution')
ax[1].legend()
ax[1].set_xlim(-4, 4)
ax[1].set_title('Histogram of Standardized Residuals')
ax[1].set_xlabel('Residuals')
ax[1].set_ylabel('Density')
plt.tight_layout()
plt.show()

Walidacja założeń modelu

Autokorelacja residuów - ACF

from statsmodels.graphics.tsaplots import plot_acf

fig, ax = plt.subplots(1, 1, figsize=(14, 4))
plot_acf(model_fit.std_resid, lags=40, ax=ax, alpha=0.05)
ax.set_title('ACF of Standardized Residuals')
plt.show()

Walidacja założeń modelu

Test Ljung-Boxa

Test Ljung-Boxa sprawdza, czy reszty modelu są niezależne. Hipotezy testowe są następujące: \[ H_0: \text{reszty są niezależne} \\ H_1: \text{reszty są skorelowane}. \] Robimy wykres \(p\)-value dla różnych opóźnień \(h\).

from statsmodels.stats.diagnostic import acorr_ljungbox


ljung_box = acorr_ljungbox(model_fit.std_resid, return_df=True)
ljung_box['lb_pvalue'].plot(title='Ljung-Box Test p-values', figsize=(14, 4), style='.')
plt.axhline(y=0.05, color='r', linestyle='--')
plt.xlabel('Lag')
plt.ylabel('p-value')
plt.show()

Miary dopasowania modelu

Akaike Information Criterion (AIC)

\[ AIC = -2 \cdot \log(L) + 2 \cdot k \] gdzie \(L\) to funkcja wiarygodności, a \(k\) to liczba parametrów w modelu. Im mniejsza wartość AIC, tym lepsze dopasowanie modelu.

  • Preferuje modele o mniejszej liczbie parametrów, ale nie karze ich zbytnio za złożoność.
  • Jeśli zależy nam na lepszych zdolnościach predykcyjnych, to lepiej użyć AIC.

Bayesian Information Criterion (BIC)

\[ BIC = -2 \cdot \log(L) + k \cdot \log(n) \] gdzie \(L\) to funkcja wiarygodności, \(k\) to liczba parametrów w modelu, a \(n\) to liczba obserwacji. Podobnie jak AIC, im mniejsza wartość BIC, tym lepsze dopasowanie modelu.

  • Preferuje modele o mniejszej liczbie parametrów, ale karze je bardziej za złożoność.
  • Jeśli zależy nam na lepszym dopasowaniu modelu do danych, to lepiej użyć BIC.

Backtesting modelu GARCH(1,1)

Miary

Mean Absolute Error (MAE)

\[ MAE = \frac{1}{n} \sum_{t=1}^{n} |y_t - \hat{y_t}| \]

Root Mean Squared Error (RMSE)

\[ RMSE = \sqrt{\frac{1}{n} \sum_{t=1}^{n} (y_t - \hat{y_t})^2} \]

  • Wykonujemy prognozę o jeden krok do przodu, a następnie porównujemy prognozowane wartości z rzeczywistymi wartościami. Wykonujemy prognozę na podstawie modelu GARCH(1,1) i obliczamy MAE i RMSE.
from sklearn.metrics import mean_absolute_error, mean_squared_error

# Wykonujemy prognozę na 1 krok do przodu
forecast = model_fit.forecast(horizon=1)
predicted = forecast.variance.values[-1, :][0]
actual = googl.iloc[-1]**2
mae = mean_absolute_error([actual], [predicted])
rmse = np.sqrt(mean_squared_error([actual], [predicted]))
print(f'MAE: {mae}')
print(f'RMSE: {rmse}')
MAE: 0.008448250442643305
RMSE: 0.008448250442643305

ARMA-GARCH w pakiecie R

#| echo: FALSE
#| include: FALSE

library(xts)
library(zoo)
library(rugarch)
library(ggplot2)
library(dplyr)
#| echo: FALSE
#| include: FALSE

googl <- read.csv('data/googl.us.txt', header = TRUE, sep = ",")
googl_xts <- xts(googl[, 5], order.by = as.Date(googl[, 1]))
log_returns <- diff(log(googl_xts))
log_returns <- na.omit(log_returns)
ggplot(log_returns, aes(x = index(log_returns), y = log_returns)) +
  geom_line() +
  labs(title = "Log Returns of Google Stock", x = "Date", y = "Log Returns") +
  theme_minimal()
#| echo: TRUE
spec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
              variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
              distribution.model = "norm")
fit <- ugarchfit(spec, log_returns)
print(fit)
ggplot(data.frame(fitted = fit@fit$fitted.values, time = index(log_returns)), aes(x = time, y = fitted)) +
  geom_line() +
  labs(title = "Fitted GARCH Model", x = "Date", y = "Fitted Values") +
  theme_minimal()
ggplot(data.frame(residuals = fit@fit$residuals, time = index(log_returns)), aes(x = time, y = residuals)) +
  geom_line() +
  labs(title = "Residuals of GARCH Model", x = "Date", y = "Residuals") +
  theme_minimal()